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ABSTRACT 

This Thesis develops an efficient and effective method for designing and 
analyzing the performance of various integrated optical waveguide structures using the 
beam propagation method of analysis. Modifications in the physical layout of an optical 
device through changes in coupling connection design, splitting angles and waveguide 
dimensions may have significant effects on device performance. The beam propagation 
method is initially developed for a symmetric Mach-Zehnder interferometer for baseline 
validation of the accuracy and applicability of the propagation scheme. A major 
validation is achieved through modeling an asymmetric device designed and built by the 
Naval Research nen The validated simulation model is used to analyze the 
performance and design characteristics of complex parallel configurations of 
interferometers. The beam propagation method allows quantitative analysis of the 
performance of these integrated optical devices. The propagation model developed 
implements a new global propagator scheme that substantially reduces computational 
requirements and introduces a design methodology that ensures compatibility between 
the discrete implementation and the physical structure. Also identified are areas in which 
continued research can provide a complete modeling system that may be implemented as 


a stand-alone design and analysis tool. 
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L INTRODUCTION 


A. BACKGROUND 

A number of computational schemes exist which can be used to numerically 
model integrated optical waveguide structures. These schemes include the second-order 
finite differencing method, the split operator method, the short iterative Lanczos 
propagator and the Chebyshev scheme. These schemes have been utilized to propagate 
solutions to the scalar Homogeneous Helmholtz equation and the time dependent 
Schrédinger equation as well [1]. The applicability of the propagation scheme to the 
physical structure, numerical efficiency, stability of the solution and computational 
constraints must all be considered when choosing a modeling scheme. 

The beam propagation method (BPM) is a split operator method to solve the 
scalar Homogeneous Helmholtz equation. The BPM has been shown to be an effective, 
efficient tool for generating required numerical solutions for waveguide structures [1]. 
The BPM permits the rapid analysis of design parameters and overall system performance 
in an extremely short period of time. 

B. PRINCIPLE CONTRIBUTIONS 

Although the BPM has been previously utilized for optical analysis [2, 3], the 
time constraints and required flexibility for analyzing relatively long complex integrated 
optical devices posed serious Sebistns, This thesis develops an optimal BPM algorithm 
for integrated optic circuits. This thesis develops a methodology for ensuring that the 


discrete parameters used to implement a physical device are matched to the device 


parameters by utilizing a subtle aspect of sampling theory that forces equalization. A 
major contribution of this thesis is the development of the global propagator scheme. 
The global propagator developed in this thesis has made the BPM an important tool for 
' designing and analyzing the long structures required for highly sensitive voltage detection 
devices such the Mach-Zehnder Interferometer (MZI). 

This thesis attempts to prove the applicability of the BPM through modeling an 
existing device and comparing the theoretical BPM results to actual data measured in the 
_ laboratory. The resulting confirmation of the validity of the BPM demonstrates the 
potential applicability of the method. . The validated BPM is utilized to analyze more 
complex parallel MZI structures. The performance of parallel structures is evaluated as a 
function of separation distance between the individual interferometers, and the effects of 
radiation mode coupling are quantified. 

C. THESIS ORGANIZATION 

This thesis begins in Chapter II with the mathematical development of the BPM 
algorithm. The BPM algorithm is a numerical solution to the scalar Helmholtz equation. 
The equation is broken down into an alternate series of propagation and lens terms that 
allow the optical field to propagate through a guiding structure. The constraints and 
limitations of the BPM are analyzed and a complete set of algorithm design equations are 
developed. This development provides the basis for the global propagator scheme that is 


employed in this thesis. 


tN 


In Chapter [II the first test of the BPM is the simulation of the well known 
symmetric MZI. Since the single symmetric MZI has been previously analyzed the 
expected output is available for comparison. Although the equations for the output of the 
single MZI do not include radiation losses, this first step is vital in developing a reliable 
design tool. This provides the basis for validation of the global propagator, and 
demonstrates the methodology for implementing the discrete matched device parameters. 

The modeling of the single asymmetric MZI is a major validation phase of the 
code development performed in Chapter [V. The BPM representation of the asymmetric 
interferometer is developed and the BPM analysis is compared to the physical data 
measured in the laboratory. The BPM is shown to be a valid integrated optical design 
tool that effectively implements radiation losses due to device characteristics. 

The complex analysis of a parallel configuration of Mach-Zehnder Interferometers 
is performed in Chapter V. Parallel configurations of guided wave MZI are important in 
signal processing architecture that utilize voltage modulators. The BPM is used to 
analyze the effects of radiation mode coupling between adjacent interferometers and 
radiation losses due to branching angles and power dividers. 

Finally, Chapter VI addresses the potential use of the BPM for design of other 
integrated optical devices. The use of integrated optics is a growing area of signal 
processing and the use of optical computers is rapidly becoming a reality. Potential 
improvements of the BPM algorithm to include other optical effects and unusual design 


parameters are addressed. 


Il. DEVELOPMENT OF AN OPTIMAL BPM ALGORITHM 


A. BEAM PROPAGATION METHOD 


The Beam Propagation Method (BPM) is a numerical solution to the wave 
equation that effectively models an optical structure as a series of infinitely thin lenses 
separated by an incremental axial distance Az. 

1. Fresnel Approximation to the Wave Equation 

If we assume that light propagating in an optical waveguide is monochromatic, 
then we can apply the homogeneous Helmholtz Equation 


2 
ey: o 
VE +7170, x, )E = 0, (1) 
where the refractive index n?(@, x,y) is assumed to only have dependence on the x and y 
coordinates [4]. If we assume a forward propagating field in the + z direction, then we 


can express Eo, x, y,Z) as 7 
E(@, x,y,z) = Eexp(-jk2), 


(2) 
where 
k = now/c, G) 
and n, is the refractive index of the substrate. 
By substituting Equation (2) into Equation (1) we now have 
: 2 
V?Eexp (~jkz)+ er Gs y)E exp (-j kz) = 0. (4) 
Since 
oC? o2 eo 
Ves ste tes, 
ox? ay? az? () 
we can evaluate the partial derivative with respect to z separately, so that 
Pe 6 [ 8 F | 
Aen jkz)= 25) a Eexp(-jkz) (6) 


= 2 exp, jkyZE jkBexp(-jkz) | (7) 


2 
= exp ( jk2) SSE jk exp ( jkz) 2B -kexp( jk z)E - jk exp ( jky Ze (8) 


= exp (-jk2y 5 oR 2 jk exp (-j kz) = oR k?exp ~jkz)E. 9 
< ikz) (9) 


By substituting Equation (9) into Equation (4) we now have 
a é 


n? 
a5 it is 2jk aE - WE+Vin+k) GE = 0, (10) 
where 
cnn Pace 
Viz Be * By? (11) 
Now, rearranging the terms of Equation (10) into a more recognizable form we have 
2 pir, Ope? (a 
a ar es mot E. (12) 


If the field is slowly varying in the + z direction then we can neglect the first term on the 
left in Equation (12). This gives the paraxial or Fresnel form of the wave equation 


2 
2jk SE! = vip’ ie 82-1} (13) 
9 


or 
an, J 2.42 m1] hex 14 
= E rd jviek e 1| |E’ =0, (14) 


where E’ represents solutions to the Fresnel approximation to the wave equation. Since 


Equation (14) represents a standard ordinary differential equation, the solution is 


E’(x, y, Az) = oxp| -L22[ vi « oe 1}) rex. (15) 


where E/(x, y,0) is the initial condition of the field at z = 0. Separating the operation of 


the exponent in Equation (15) provides the basis for the Split Operator Method (SOM) or 


BPM [2], so that we now have 


( jAz—o\ f faz afn2 : 
E’(x, y, Az) = exp) -=— Vii jexp| -=— k?| == - 1] | E(x, y, 0). (16) 
\"2k err: n2 ( 


As developed in [2], Equation (16) can be further split so that 


jAzk 
E’(x, y, Az) = exp(—12v3 Je xp | iae ea )) exp( LZ vt }B’(xy.0) 
+O(Az)’, 


(17) 
where the three exponential factors represent a half step of propagation, a phase or lens 
term, and a second half step of propagation respectively. Since the V? term in the first 
exponential term of Equation (16) does not commute with the x and y dependence of the 
second exponential term, an error term OfAz)’ is introduced in Equation (17). The error 
term is not evaluated at this point since the split operation is not actually performed in the 
BPM implementation. In structures that require large axial distances, adjacent 
propagation steps are normally combined to reduce computation time, as will be shown in 
a later section. However, Equation (17) provides an excellent method of visualizing the 
propagation process by replacing the optical waveguide with a series of infinitely thin 
lenses [1], as shown in Figure 1. 

If the initial field is composed of a limited spectral bandwidth then the field can 


be represented by a truncated Fourier series such that 


! _ we eA 2nmx | 2nny \ | 
en a coe BF mn(2) exp Poy Ly L2 ‘he ay g®) 


where L, is the width of the computational grid in the x direction, L, is the width of the 


computational grid in the y direction and N and M are the total transverse points in the 


Figure 1. Thin lens representation of BPM algorithm. 


computational grid. If the propagation term 
j Az 
exp (- ae (19) 


operates on the Fourier series representation of Equation (18) the result is that 


ak ae J OAS 


This technique allows the field to propagate forward by calculating new Fourier 


E’ mn(Z + Az/2) = E’ mn(z)exp iz (220) ; + (2x8 | ‘ | (20) 


coefficients based on the spatial frequency propagator 
jAz((2nm)* (2=2)') 1 
vp (22 } oad (21) 
that is used in Equation (20). In effect this leads to the determination of a new Fourier 
series representation of E’ (x, y,z+Az/2 utilizing Equation (18). The new representation 


of the complex field is then multiplied by the lens term 


exp i & - 7 (22) 


of Equation (17). This in turn leads to the determination of another Fourier series 


representation of the field at E’ (x,y,z+Az/2 immediately after the lens, and the 


propagation step is applied once more. The application of Equation (20) allows the use of 
Fast Fourier Transform algorithms to implement the propagation terms in Equation (17). 
This concept allows for a much easier algebraic implementation of the propagation term 
in the spatial frequency domain, while the phase or lens term is performed with simple 
multiplication in the spatial domain. Thus, the BPM consists of a set of Fourier 
transforms interspersed with complex multiplications in the spatial domain in an iterative 
algorithm that advances the solution in successive steps along the optic axis [5]. The 


Fourier transforms are carried out using the well known Danielson-Lanczos FFT 
algorithm and require transverse grids of 2™' x 2™ in the spatial and spectral domains, 


where m, and m, are determined by N and M of Equation (18) and are subsequently 
analyzed. 
2. Direct Solution to the Wave Equation 
An alternate method of developing a split operator solution to the wave equation 
without initially making the Fresnel approximation has also been utilized [2]. If Equation 
(1) is rewritten in the form 
cE + (vt + on, )) E=0, (23) 
the representation of Equation (23) can be treated as a second order ordinary differential 
equation. The solution at z= Az may be written in terms of the field at z = 0 as 
E(x, y, Az) = exp E ad v3 + 22] "Tee. y, 0). (24) 


2 
c 
The square root in the exponent of Equation (24) can be written in the form 


3 5\ 12 We 
\ 12 (v3 + (onlc)?) IG + (onic)?) + (oni)| 
a (25) 


(o2, @ 
Rae Pp (v2 ae 2) 2 
Vi (wn/c) +(on/c) 


(vi + (onlc)?) + (v3 + (onle)’) Y taniss 
(26) 


(V3 +(onley?) "” + onic) 


2 (wnley? +(onle(V3 + (nlc) i 
a 
= Q7) 


ear 1/2 
Vii +(onlc)’) + (on/c) 


(v3 if (onlc)*} iat (onic) 


(onic) (oni + (v3 # (onio)) i 
(28) 


vi 
i2 
(v2 + (onlc)?) + (nic) 


ae 12 
(v3 + (onic?) +(a@nic) 


Vi +(onle). (29) 


(v3 + (onle)’) i + (an/e) 


If the variations in n(x,y) are sufficiently small, then the » in the first right hand member 


of Equation (29) can be replaced with n, and utilizing Equation (3) we have 


2 
(vi+%n) sheet (30) 
(vise) +k 
2 
= Vik kk G1) 
(v2.40) +k 
v2 
7 +k+k( 2-1) (32) 


If we now substitute Equation (32) into Equation (24) and assume a + z propagating 
wave, we can apply Equation (2) and we have 


f \ 


vy? 
E(x, y, Az) = exp | ~j Az— <8 
\Vi+ 2) +k 


f 
+k = 1) E(x, y, 0). (33) 


At this point the resemblance to the Fresnel approximation solution given in 
Equation (15) is readily apparent. However, if V2. in the denominator of Equation (33) is 


neglected in comparison to &’, an almost exact reproduction of Equation (15) is 
recovered. The representation of Equation (33) can now be written as 
Ex, y, Az) = exp ( ity? | exp| ~jkAe( 2 = 1) cs y, 0). (34) 
Using the same split operator technique that was used in developing Equation (17) 
E(x, y, Az) = exp e ey exp [i es _ 1) lex xp (- ity? E(, y, 0) 

+ O(Az)?, (35) 
where the same error term O(Az)’ used in Equation (17) is used here due to the 
commutation error of the transverse Laplacian. This expression is now seen to be almost 
identical to the split operator developed in Equation (17). The impact of the change in 
the lens parameter on optical structures is evaluated in subsequent sections. The only 
difference between the split operators of Equations (17) and (35) is a change in the lens 


term from 


exp ieee (m2 - i (36) 


ll 


to 
exp [race 2 - 1) | (37) 

However, it must be noted that approximations were made in deriving both 
- equations (30) and (33). Subsequently, the Fresnel approximation leading to Equation 
(17) is used as the primary algorithm in this analysis. 

3. Effective Index Method 

One of the main objectives of this thesis was to develop and analyze the 
_ applicability of a fast, efficient BPM algorithm that is adaptable to varying structures. 
The first step in this process is to reduce the three-dimensional optical structure that is to 


be analyzed to a two-dimensional structure. Analyzing Equation (18) it is readily 
apparent that a two-dimensional (44x N) point FFT must be carried out for each 


propagation step over the axial distance Az/2. Assuming an effective step index change 
in the guiding structure allows the overall device to be reduced to one transverse 


component in x, so that the refractive index can be written as 


n(x) =n + 8n- rect(% =%e) F (38) 


where 5n = ng — np is the index difference between the waveguide and the substrate and w 


is the width of the waveguide [3]. This obviously results in a reduction of M FFT's for 
each propagation step over a potentially long structure. 


Implementing a one-dimensional cross-section results in a reduction of Equation 


(17) to 
E’ (x, Az) exp 122 exp — (= 1) Jexp (LES Jee) 
+O(Az)°, (39) 
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and a reduction of Equation (18) to 


N/2 
) = t [ 27mx 
BGuz)= % Bh @) exp| (22H) | (40) 


where L is the width of the computational grid in the transverse x direction. The number 
of gridpoints N must be determined based on physical constraints of the optical system. 
The number of grid points and therefore the transverse sampling interval Ax a 
determined for each structure and are developed in their respective sections. 
The same concept that was used to propagate the three-dimensional solution is 
used on the two-dimensional solution. If the propagation term 
XP | —"Ty a? (41) 
of Equation (37) operates on the Fourier series representation of Equation (38) the result 
is that 
El m(z-+Az/2) = E! m(z)exp | (2am) ‘| (42) 
Once again, this technique allows the field to propagate forward by calculating new 
Fourier coefficients based on a spatial frequency propagator 
exp iS (22)'] (43) 
One of the major contributions of this algorithm was developed at this point. 
Noting that the spatial frequency propagator depends upon physical parameters of the 
system, but not on the lens structure, the propagator can be determined at the initial stage 
and stored as a propagator array. If the propagator is predetermined it reduces the 
number of FFT's required in the two-dimensional structure by the number of 


propagation steps in the overall structure. The total number of steps in an optical 
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structure is typically on the order of several thousand steps, as is shown in the analysis 
performed in subsequent sections. 
4, Algorithm Implementation 

a. Absorption Window 

The optical devices analyzed in this thesis are modeled using a step index 
profile, as mentioned previously. The step index profile is potentially problematic due to 
the discontinuities in the lens structure [3]. These discontinuities in the index profile 
excite spurious high spatial frequencies (i-e., radiation modes, and therefore energy 
leakage) when the BPM analysis is carried out [5]. The optical devices analyzed are 
primarily interferometers, but also include variations of Y-junction power dividers 
(YPD's). Interferometers inherently incorporate various branches in the guiding 
structures, as will be shown. These junctions and branches contribute significantly to the 
overall radiation losses of the device, and are therefore a major source of energy leakage. 

Due to the periodic nature of the Fourier transform and the finite structure 
of the system, as these radiation modes propagate to the edge of the computational 
window in the transverse plane, they are folded back to the opposite edge of the window 
in subsequent propagation steps. This energy is seen as high frequency noise and may 
cause high frequency numerical instabilities [5] depending on the physical system being 
modeled. To avoid this potential problem, the radiation modes are absorbed at the edge 
of the window. This is seen as a valid approach since the radiation loss itself is calculated 


by measuring the contained power in the guiding structure. 
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The absorption window is implemented using a gradual field absorber near 


the window boundaries given by [2. 6] 


1, Ix] < |xal 
absorb(x) = 4 1/20. + cos[m(v—xa)/ta—x5)]), — [xal <Ixl |x], (44) 
0, lxal < [xl <[x-| 


where x, is the coordinate of the grid boundary, x, represents the inner edge of the 
absorber, and x, is the outer edge. The parameters of the absorbing window must be 
chosen for each device independently to ensure the field is absorbed over a sufficiently 
large region while ensuring that the absorber itself does not interfere with the guided 
modes of the device. 

Although the absorption window is tailored for each device, as mentioned 
earlier, a generic absorber utilizing Equation (44) is shown in Figure 2 for illustrative 
purposes. In the example in Figure 2, the parameters x, = 200, x, = 255, and x, = 256 are 
used. 

b. Transverse Sampling Interval 

An additional problem introduced by discretizing an optical circuit is 
discovered in determining the transverse and axial sampling intervals. This problem is of 
course compounded when attempting to implement a discontinuity such as a step index 
profile. Five factors must be considered [3, 5] when determining the transverse sampling 
interval Ax: (1) the size of the FFT window, (2) the range of the field's angular spectrum, 


(3) computer limitations, (4) the tolerable trapezoidal distortion of the step index caused 
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Figure 2. Typical absorption window. 
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by the finite grid, and (5) the tolerable index boundaries location uncertainty, also due to 
the discrete nature of the grid. 

The first parameter determined is the window size W, which represents the 
-physical dimension of the system being modeled. The window size must include the 
entire optical structure and must allow sufficient distance for decay of the evanescent 
fields and the absorption window. Since the window is dependent upon the specific 
system, W and therefore Av must be determined separately for each device. 

The propagation of the light in the waveguide is considered to be paraxial 


since the direction of propagation is predominantly forward. The physical spectral range 
is chosen to be that which corresponds to propagation within |6| < 7/4 of the optical (z) 


axis. Utilizing this limitation on the spectral range, the highest angular frequency to be 
represented by the grid is given by [3, 5] 

K, =k sin (7/4), (45) 
where k is given by Equation (3). However, the highest angular frequency that can be 


represented by an FFT with N points is given by 
kEFT = (Qn N} 


IW (46) 
Due to the limitation of the FFT 
‘ 2nN 
Kx Sky" = (228) (47) 
Therefore the number of gridpoints V should satisfy 
Nz ato sin (w/a), (48) 
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For a given structure of window size W, the transverse sampling interval is readily found 
to be 
ax=(#), (49) 

and as previously noted, must be determined for each device modeled. 

c. Axial Sampling Interval 

AS previously mentioned the split operator development includes an error 
due to the commutation error and conditions on the maximum axial sampling interval 
must be determined. For the BPM to be applicable, the following four conditions must 


be satisfied [3, 5, 7]: 


(oft) +(g)' +2) REP) <ahtea(Q). 
vat ('et wey}. on 
Az << [125° (pts), (52) 

Az << 6k(p +s)”, (53) 


where 5n is given in Equation (38) and s and p are the highest transverse spatial 
frequency components of the band-limited index profile and the propagating field, 
respectively. The Fourier transform of the step index profile defined by equation (38) is 
given as 


3{8n rect(x/w)} = wsine(wf,), (54) 
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where f, is the angular frequency component. Since the majority of the sinc function's 
energy is contained within the first three lobes of the function (greater than 97%), p is 
taken to be the third zero crossing of the sinc function, so that 
w-fxo4, =3 (55) 
or 
p= Win) _ 3, (56) 
Since the field distribution is much smoother than the index profile (due to 
the step index), we also take the value of p to be a more stringent requirement for s, so 
that 
S=p (57) 
is a valid bound on s. Substituting the values of p and s into the applicability conditions 
of Equations (50) through (53), Equation (53) is found to be the most stringent and is 
subsequently used for determining the axial sampling distance Az. 
Now that a complete set of equations have been developed for the BPM 


algorithm and the expected limitations have been analyzed for general optical systems, 


the next step is to develop a practical implementation for the modeling system. 
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I. BPM IMPLEMENTATION OF A SINGLE SYMMETRIC MACH-ZEHNDER 
INTERFEROMETER 


A. ARCHITECTURE OF THE MACH-ZEHNDER INTERFEROMETER 

The architecture for the first system modeled is shown schematically in Figure 3 
and consists of a single MZI with equal path lengths (symmetric) and Y-power dividers 
and combiners. The corresponding parameters are shown in Table 1. The separation 
distance d, between the branches of the interferometer is varied by changing the length 
of £2 and by changing the branching angle a. The index of refraction for the step index 


waveguide is given by Equation (38), where n, = 2.2 is the LiNbO, substrate index and 
the static index difference between the waveguide and substrate 6” = 5.1x10°. The input 
waveguide has a width of w, = 3.0 um and the maximum guide width in the YPD is given 
by w, = 5.61 um. The branching angle a is initially set at 1°, but is subsequently varied 


in the analysis in order to quantify the impact of branching angles on radiation losses. 
1. Linear Electrooptic Effect 
The linear electrooptic (Pockels) effect provides a change in the refractive index 


proportional to the applied electric field, which for LiNbO, corresponds to [8] 


wr 
an=-[ Er) 5 (58) 


where the electrooptic tensor element r,; = 30.8 x 10°? (m/V). A typical electrode 


configuration [8] for an integrated optic device is shown in Figure 4. If a voltage V is 
applied to the electrode shown in Figure 4, the resultant electric field through the 


waveguide has an approximate magnitude of 


ial =(4), (59) 
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Figure 3. Schematic diagram of a single symmetric Mach-Zehnder Interferometer. 


22 


Electrodes 


Figure 4. Typical electrode configuration. 
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where G is the physical gap or separation between the electrodes, as shown in Figure 4. 
However, neither the applied field nor the optical field is uniform. The effective 
electrooptically induced index change within a cross section of the waveguide is [8] 
Ang(v) = — ng 33 (X) Y, (60) 
where I is the overlap integral between the applied electric field and the optical field. [ 
is defined as 


P=S [JE acc(x,y)|Eoptia(x,y)|” dx dy (61) 


Table 1. Schematic parameter values. 
For an electrode gap/mode width ratio of 1, the resultant overlap integral T = 0.5 [8]. The 
electrooptic effect can now be added to the lens equation for the region where electrodes 
are to be utilized, such as the LS region of Figure 3. The total phase shift induced in the 


waveguide over the electrode length L5 is 


ABLs= ang ts (¥)r{Es). (62) 
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If a push-pull configuration such as Figure 3 is utilized. the output intensity is 


given by [8] 


4 
a 3{ Ng t33 Ls oo 
Po =cos [Aaetslar ¥]. (63) 
This gives an excellent theoretical prediction for the symmetric interferometer with which 


the BPM results can be compared. The theoretical output intensity of Equation (63) gives 


a prediction of the periodic output intensity that has a first null at 


_ AG 
: 2L5 Tr33 nz (64) 
which for this device is calculated as 
900 x 10° 10-6 


Vaz = = =8.17 
2(1000 x 107®)(0.5)(30.8 x 1071?)(2.2)° 
This would produce a voltage folding period of approximately 16.34 V, which can be 


used to verify the accuracy of the BPM. However, Equation (63) does not take into 
account radiation losses due to branching angles or power dividers. 

2. Structural Parameters 

In developing structural parameters, and therefore BPM implementation 
parameters, the guidance or mode confinement characteristics of the waveguide must be 
considered. The guidance strength in a three-layer waveguide is determined by its 
normalized frequency V,,, defined by [5] 

Vi= =t8 (2 ny 6n). (66) 

A higher V,, means stronger guidance and therefore better mode confinement. For a given 
V,,, a waveguide can support M guided modes, where 


(M-1)a< Va< Ma. (67) 
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For this device A = 0.9 um was used. Since n, and 6n have already been determined, the 
choice of w, determines V,. In this analysis strong guidance with single mode operation 
was desired, therefore w, = 3.0 4m was chosen. This results in a normalized frequency V, 
= 3.137, which provides single mode operation. The next step in implementation is to 
determine the overall grid size. The grid size is needed in Equation (40) for the spatial 
frequency propagator and is also utilized in Equation (46) when determining the 
transverse sampling interval. In the initial analysis, the arm separation distance d, shown 
in Figure 4 is forced to be 20 um by choosing the appropriate length for L4, with a fixed 
branching angle a. If an absorption window of 7 1m is implemented on each side of the 
grid and a radiation mode decay distance w,,.., is chosen as 10 um, we can determine an 
appropriate grid size as 
Ws dgt+2w1+2 Weecay +2 Wabsors = 60 Lm. (68) 

In this analysis the MZI structural parameters will be manipulated so an additional buffer 
of 20 um is proposed to be added to Equation (68), resulting in W = 80 um. 

3. Sampling Interval 

At this point, all of the information required to determine the transverse sampling 


interval Ax as given by Equations (48) and (49) is available. Since we know from 


Equation (48) that 
Ne a M sin (n/4), (69) 
we therefore have 
2(2.2)(80 um) (_1_) _ (70) 
Ne o00nm = (ia) 27 
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Since the algorithm is to be implemented using FFT's, as previously mentioned, this 
would require NV = 512 for this circuit. Noting that W = 80 um provides much more 
radiation mode decay distance than is required by Equation (68), if W = 69.81 um is 
‘chosen, an additional 9.81 um additional distance over the Equation (68) requirement is 
still provided, the significance of which will be shown directly. Applying W = 69.81 ym 


to Equation (69) we now have 


N2 241, 
(71) 


so N = 256 is chosen. The significance of the relationship between N and W is shown by 
utilizing Equation (49), where 


69.81 
= (2-81 um | = 0.2727 um. (72) 


Ax 
Since the guide width must be implemented in discrete steps and has been chosen as 3.0 
jum, we can determine the actual discrete guide width w, by finding the number of grid 
steps 


wae Ce) | Ax = 2.996 um, (73) 


where | | denotes an integer value. This demonstrates the importance of choosing the 


correct relationship between W, N, Ax, w, and w,. As an example, if the grid size of W= 
60 um developed in Equation (68) had been used, then the result would be Ax = 0.234 um 
and using Equation (73), wy = 2.8 um. This would reduce the guidance strength of the 
model structure. Therefore the parameters chosen for this analysis are: N = 256, W = 


69.81 wm, and Ax = 0.2727 um. 
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Since the guide width has been chosen to be 3.0 tm, and the discrete guide width 
w, has been designed to be approximately 3.0 tum as well, the requirements for the axial 


sampling interval can now be determined. Solving Equation (56) we find 


_{ 3m) oom 4 
p=| ~ ym? ( ) 


3 um 
and as noted in Equation (57), p = s. As previously noted, the most stringent criterion is 
given by Equation (53) so that 
Az << 2.4um (75) 

In order to reduce computation time Az = 2.4 um is used in these simulations, which 
provides valid results when compared to the theoretically expected values. All of the 
structural parameters for the BPM implementation have now been developed. The last 
factor to be determined is the mode distribution of the optical input field. 
B. SYMMETRIC EIGENFUNCTION INPUT 

The input eigenfunction u,(x, z = 0) to the interferometer is computed analytically 


using an algorithm implemented in a C™ program. The normalized field distribution in 


the input waveguide is given by [3, 5] 


_ my _ | 98{Ki0 x), Ix] Swi/2 (76) 
Malt, z= O)= cos(k1o w1/2)exp [—k20(Ixl — wi/2)],. xl 2 wi/2 
where 
sen 2 1/2 
Kw=|( ne 8 | (77) 
and 


5 a2 
k= | 03-( uit | . (78) 


The fundamental mode eigenvalue B, is extracted from the transcendental eigenvalue 
equation 
=T (79) 
For the 3.0 um waveguide used in this analysis B, = 15381674 cm". The input field that 
is to be launched in the BPM algorithm is calculated numerically. Figure 5 depicts the 
symmetrical mode eigenfunction that is launched into the guiding structure. 
C. BPM IMPLEMENTATION 

1. Program Structure and Implementation 

Now that all of the parameters of the BPM algorithm have been determined, the 
structure to be simulated has been chosen and the input field has been selected, the final 
step is to implement and demonstrate the BPM. A detailed flow diagram of the BPM 
algorithm is shown in Figure 6. As mentioned previously, one of the key contributions of 
this thesis was the reduction in computational time produced by implementing a global 
propagator. Since the propagator is applied in the spatial frequency domain, an analysis 
of the optical field in the frequency domain provides invaluable insight into the operation. 
The optical field is stored in an array during each propagation step along the optical axis 
and the FFT is applied to the optical field. The Fourier transform of the optical field 
demonstrates the symmetry that is intuitively expected and is shown in Figure 7. Another 
major reduction in computation time is achieved by noting the effect of this symmetry of 
the Fourier transform. If the propagator array is prefolded then it can be directly 


multiplied by the optical field in the frequency domain, eliminating the need for bit 


29 


Normalized symmetrical field 


1 i 
| A | 
| i 
if 
0.8- [| 
; s | 
poy : 
| | | 
0.6. [| 
} P| 
| 
0.4. i | 
! ! i 
or 
i 1 
i ‘ 
0.2. a 
i if i 
: \ 
! / \ 
/ \ ; 
0. pital aces af et n _ Sis 


0 50 100 150 200 


x in gridpoints 


Figure 5. Normalized symmetric eigenfunction input. 


30 


Absorb(x) 


Structure 
? 


Figure 6. Detailed BPM flow diagram. 
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Figure 7. Fourier transform of the symmetric eigenfunction input. 


32 


reversal during the propagation loop. The normal and prefolded complex propagators are 
shown in Figure 8(a) and 8(b) respectively. 

The beam propagation method of analysis provides a convenient method of 
viewing the launched eigenfunction as it propagates down the optical structure. The 
BPM implementation of the guide structure shown in Figure 3 is detailed in Figure 9(a). 
To demonstrate the radiation modes being studied, Figure 9(b) details the magnitude of 
the optical field [E'(x,z)| as it propagates down the interferometer. The spacing between 
the arms of the MZI is d, = 20 um. From these figures the periodic radiation mode 
coupling between the arms of the interferometer is readily apparent. Figure 10 shows a 
cross-sectional view of the center of the interferometer superimposed upon a scaled 
model of the step index profile (for reference), and shows the mode distribution within 
each waveguide with 0 volts applied. A detailed view of the Y-power divider designed 
for this interferometer is shown in Figure 11. The BPM implementation of the index 
structure is shown in Figure 11(a) and the BPM calculated optical field distribution is 
shown propagating through the structure in Figure 11(b). A key insight into the 
mechanics of the MZI can be obtained by viewing the output characteristics. Figure 12 is 
a detailed view of the output power combiner for this device. The output characteristics 
are demonstrated in Figures 13(a) and 13(b), with 0 volts and 8.17 volts applied to the 
electrodes receptively. The theoretical null predicted by Equation (65) is verified through 
the demonstration in Figure 13(b). The impact of the modulation voltage is clearly 


demonstrated by viewing an instantaneous cross-section of the device. A cross-sectional 
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Figure 8. Global propagator: (a) normal; (b) prefolded. 
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(b) 


Figure 9. Single symmetric MZI: (a) step index profile; (b) BPM analysis. 
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Figure 10. Cross-sectional view at the center of LS with V = 0. 


36 


(b) 
Figure 11. Input YPD design analysis: (a) step index profile; (b) BPM analysis. 
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Figure 12. Output power combiner design: step index profile. 
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Figure 13. Output power combiner BPM analysis: (a) V = 0 volts; (b) V = 8.17 volts. 
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view of the output section 27 is shown in Figure 14(a) and 14(b) for 0 volts and 8.17 
volts applied. This view clearly shows the desired interference effect introduced by the 
modulating voltage 
The BPM propagation algorithm is implemented in a C™ program. The 
implementation allows branching angles, structural lengths and electrode voltages to = 
input through external script files. This enables rapid assessment of the impact of 
changing design parameters on the overall device performance without recompiling. This 
gives the potential for a stand alone, platform independent design package. A typical 
flow chart of a BPM analysis showing the usual sequence of implementation and analysis 
is shown in Figure 15 As shown in Figure 15, once the eigenfunction is calculated for a 
specific structure it is stored for repeated analysis. A command line script file is utilized 
to provide structural or electrode voltage parameters so that detailed analysis of system 
performance can be performed without continual user input. A set of script files used for 
manipulating voltage and structure parameters for the devices in this thesis are shown in 
the Appendix. 
2. Validity and Applicability of the BPM 
The theoretical output intensity of the interferometer being modeled is given by 
Equation (63). Varying the electrode voltage from 0 to 20 volts and comparing the BPM 
output intensity to the theoretical intensity of Equation (63) gives a good analysis of the 
validity of the BPM. The BPM calculated output intensity is shown versus the theoretical 


output intensity in Figure 16. The 16.34 volt folding voltage predicted by Equation (65) 
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Figure 14. Cross-sectional view of the output: (a) V 


0; (b) V = 8.17. 
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Figure 15. Flow diagram of BPM analysis sequence. 
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Figure 16. Theoretical versus BPM calculated output intensity. 
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is in fact closely tracked by the BPM analysis. The results of Figure 16 imply that the 
BPM provides an accurate prediction of interferometer performance. The alternate lens 
equation developed in Equation (35) was also used to analyze the symmetric MZI. As 
’ demonstrated in Figure 17, the alternate lens equation has little effect on the analysis. 
Design parameters of an MZI can have dramatic effects on optical device system 
performance, and can be easily modeled using the BPM. Varying the length of section 
L4 in Figure 3, while holding the branching angles constant, changes the arm separation 
_ distance d,. The impact on radiation mode losses and coupling can be significant. The 
radiation mode coupling with arm spacing d, = 5 um is shown in Figure 18, while the 
results shown in Figure 19 demonstrate the effect of increasing the spacing to d, = 10 pm. 
It would seem that by increasing arm separation radiation mode coupling would be 
eliminated. However, a major problem in system design is performance optimization. 
As arm spacing is increased, so is the overall device length. The device length has a 
direct impact on loss, which is normally to be minimized. The overall device loss 


through the interferometer is calculated as 


( Pout | 
a 


Loss = 10 log (80) 
The device loss for arm spacing d, = 5 xm, 10 tum, and 20 um was calculated as 1.6, 0.49, 
and 0.72 dB respectively. However, by reducing the splitting angle a of the device the 
loss can be reduced. The loss with a = 0.6 degrees was calculated for arm spacing d, = 5 


uum, 10 pm, and 20 um as 0.7, 0.1, and 0.2 dB respectively. The effects of reducing the 


splitting angle on loss are clearly evident and the reduction in radiation mode coupling 
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Figure 17. Alternate lens equation analysis of output intensity. 
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(b) 
Figure 18. MZI with spacing d, = 5 um: (a) step index profile; (6) BPM analysis. 
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(b) 
Figure 19. MZI with spacing d= 10 um: (a) step index profile; (b) BPM analysis. 
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can easily be seen utilizing the BPM. An analysis of the reduction in radiation mode 
coupling is shown in Figures 20, 21, and 22 for arm spacing d, = 5 um, 10 um, and 20 
tum respectively. The most notable difference is achieved for d, = 5 tum; comparing 
Figure 18 to Figure 20 it is apparent that increasing branching angles has a dramatic 
effect on device loss due to the lack of energy containment at increased branching angles. 
The BPM is an excellent tool for minimizing these types of loss prior to fabrication. 

The BPM has now been demonstrated as a viable tool for optical system design 
and has been validated through comparison of calculated results to theoretical predictions. 
However, as was noted previously, the theoretical prediction of Equation (63) does not 
include radiation losses and thus gives no indication of expected device loss, such as 


Equation (80). 
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Figure 20. MZI with spacing d, = Sm: (a) step index profile; (b) BPM analysis. 
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Figure 21. MZI with spacing d, = 10 ym: (a) step index profile; (b) BPM analysis. 
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(b) 
Figure 22. MZ] with spacing d, = 20 um: (a) step index profile; (b) BPM analysis. 
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IV. SINGLE ASYMMETRIC MACH-ZEHNDER INTERFEROMETER 


A. DEVICE PARAMETERS 

A true test of the validity of the BPM is achieved through modeling a physical 
device that can be tested in the laboratory. An asymmetric Mach-Zehnder Interferometer 
fabricated by the Naval Research Laboratory [9] is used for this purpose. A schematic 
representation of the physical device, showing electrode placement, is shown in Figure 
23. The representation used to develop the BPM implementation of the physical device is 
shown in Figure 24. The physical dimensions of the asymmetric device are given in 
Table 2. A major increase in computational time is intuitively expected from the overall 
device length of 21.96 mm compared to 1.32 mm for the symmetric device. 

The first step in implementing the asymmetric structure is to once again determine 
the grid parameters. Since in this case a reverse engineering process is utilized, some 
flexibility would seem to be removed. However, the discrete guide width w, must still 
match the physical guide width w,. Using Equation (68) the overall grid width 
requirement is determined to be 

W 2128 pm (81) 
However, using previously gained knowledge of the relationship between N, W and Ax, 
the grid size is chosen as W = 139.264 um. Since this device was designed to operate at a 
wavelength of 1300 nm, the number of grid points is determined using Equation (69) so 


that 


., 2(2.2)(139.264um) 1 


= 82 
~ 1300nm /2 ane e 
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Figure 23. Asymmetric interferometer showing electrode placement. 
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Lar m 


Figure 24. Schematic representation of an asymmetric interferometer. 
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Therefore, for the asymmetric interferometer we must choose V = 512 points. This 


results in a transverse sampling interval 


Ax = 0.272 wm. 


la 


Table 2. Device parameters. 
Since the guide width w, = 6.8 um is given [9], the discrete guide width is now forced to 
be wy = 6.8 um as well. The discrete guide width is calculated from Equation (73) as 
w= | SSu lax = 681m (84) 
This is achieved by choosing the appropriate value for W in Equation (68). Although 
Equation (68) requires a minimum value for W, extra range can be added to ensure the 


correct desired discrete guide width, while minimizing N. An axial sampling interval of 


Az= 2.4 um is used in this simulation as well. 
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B. BPM IMPLEMENTATION 
The BPM implementation of the physical structure is shown in Figure 25. The 
step index profile is shown in Figure 25(a), while Figure 25(b) shows the magnitude of 
. the optical field as it propagates through the device with 0 volts applied to the electrodes. 
The asymmetric MZI shown in Figure 23 has an inherent difference in the interferometer 
arm lengths, AL, which gives rise to an intrinsic phase bias $, given by [9] 
Oo = 2M Neg AL/A, (85) 
where Me is the mode effective index. The BPM algorithm is implemented using the 
same basic flow diagram shown in Figure 6, utilizing the known physical parameters of 

the device [9]. 

The output characteristics of the physical device were recorded and are compared 
to the BPM analysis in Figure 26. This gives an excellent validation of the BPM 
methodology. However, device loss must still be addressed. 

C. COUPLER DESIGN 

The Y-power divider (YPD) of Figure 23 is notably different than the YPD of 
Figure 11. One of the main advantages of the BPM is the flexibility of implementing 
various design concepts rapidly and analyzing the effect. A more detailed view of the 
input YPD for the asymmetric device is shown in Figures 27. The BPM implementation 
of the YPD index profile is shown in Figure 27(a) and the optical field distribution is 
shown in Figure 27(b). The optical field interaction at the output combiner is of major 


interest, and effects the overall dive characteristics. Due to the inherent phase offset of 
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Figure 26. (a) Measured output intensity; (b) BPM calculated output. 
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Figure 27. Input YPD analysis: (a) step index profile; (b) BPM analysis. 


the asymmetric device the first null occurs at V= 1.0 volts and the first maximum occurs 
at V = 3.1 volts. The BPM profile of the output power combiner index profile is shown 
in Figure 28 while the optical field distribution in the combiner is shown in Figures 29(a) 
and 29(b) for V = 1.0 volts and V = 3.1 volts, respectively. An easier interpretation of the 
effects at the output combiner is seen as a cross-section of the optical field is 
superimposed upon a scaled version of the index structure. The cross-section of L7 
shown in Figures 30(a) and 30(b) show the effects of an applied voltage of V = 1.0 volts 
and V = 3.1 volts, respectively. 
D. VALIDATION 

Utilizing BPM calculations of the output power characteristics, the overall 
device loss can be projected using Equation (80) as Loss = 3.42 dB. The actual device 
loss was reported [9] as approximately 3 dB. This is a significant validation of the BPM 
ability to properly include radiation mode losses due to device length and branching 
angles. Since the arm separation distance was substantial in the asymmetric device the 
actual radiation coupling between branches is minimal. However, many parallel 
processing architectures require multiple parallel devices in which radiation coupling may 
be even more significant than radiation losses in the substrate. 

Now that the BPM has been shown to be a valid tool for analysis of optical 


devices, more complex systems may be analyzed with confidence. 
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Figure 28. Output power combiner design: step index profile. 


62 


(b) 
Figure 29. Output power combiner BPM analysis: (a) V = 1.0 volts; (b) V = 3.1 volts. 
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Normalized E-field distribution 


Normalized E-field distribution 
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Figure 30. Cross-sectional view of L7 region: (a) V = 1.0 volts; (b) V = 3.1 volts. 
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Vv. PARALLEL CONFIGURATION OF MACH-ZEHNDER 
INTERFEROMETERS 


Parallel configurations of MZIs are important in many signal processing 
architectures. They can efficiently couple wideband rf signals into the optical domain 
making them useful for such applications as high-speed analog-to-digital converters, 
temperature sensing, and real time optical correlators [5]. In these types of architectures 
the interferometers are arranged in parallel configurations containing many successive 
bends and branches that distribute the optical power to the various sections of the circuit. 
A significant problem with the parallel integration of these devices is the limited 
efficiency due to radiation loss [5]. The BPM is a useful tool in designing and analyzing 
these types of systems. 

A. DEVICE PARAMETERS 

The initial architecture is shown schematically in Figure 31 and consists of two 
parallel interferometers that employ YPD's and combiners that have been analyzed earlier 
in this thesis. The corresponding parameter values are shown in Table 3. The separation 
distance d, between the interferometers is varied by changing the length of L2 while 
holding the branching angle constant. The individual interferometers are configured 
identically to those analyzed previously. However, the overall grid width must be 
increased to accommodate the parallel array. The window size is again calculated using 
Equation (66); W = 106 um would allow for the decay distances, absorption windows, 


separation distances d, = 20 um, and d, = 20 pm. Once again, the impact of the window 
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Figure 31. Schematic diagram of two parallel MZIs. 
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size must be evaluated in relation to the desired discrete guide width w,. In this case the 
guide width is desired to be 3.0 um, so the appropriate adjustment must be made to the 
grid size W. It is anticipated from previous evaluations that the number of grid points 
will be 512, due to the grid size, so W = 102.4 um is chosen. Using Equation (48) we 


solve for NV as 


2(2.2)(102.4um) | 
900mm 


and as expected N = 512 is chosen. The next step is to determine the transverse sampling 


N2 


= 354, (86) 


interval using Equation (49) 


(87) 
Table 3. Schematic parameter values. 
The discrete guide width w, can now be verified utilizing Equation (71) we have 
3.0um 
=| ——— = 3. 2 88 
Wg | 222 [o.aum 3.0um (88) 
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The validation of the guide width is important in ensuring the correct mode distribution 
and guidance strength. The axial sampling interval is maintained at Az = 2.4 um for this 
analysis. 
-B. BPM IMPLEMENTATION 

The BPM implementation of the parallel array index structure is shown in Figure 
32(a) and the optical field distribution is shown in Figure 32(b). In the analysis shown, 
the spacing between the interferometers d, = 10 ism and the radiation losses of the 
individual interferometers is the same as previously noted. However it is apparent from 
the optical field distribution the radiation mode coupling between the adjacent 
interferometers does exist. A clear view of the radiation mode coupling is given in the 
cross-sectional of the ZS region shown in Figure 33. The cross-sectional view shows the 
optical field distribution super-imposed upon a scaled version of the waveguide index 
structure and clearly demonstrates the existence of power coupling. 
C. RADIATION MODE EFFECTS 

The radiation modes shown in Figure 32(b) may extend to the point that 
interference with adjacent interferometers exists. The effects the radiation modes have on 
the interferometer outputs and the respective change in the mode propagation constant 
(AB) are a function of the separation distance d,. The radiation modes of one 
interferometer may also affect the radiation modes of the other interferometer. This 
interference effect produces a finite AB within the affected interferometer. The radiation 


can be eliminated by extending the separation distance d, to infinity. However this is 
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(b) 
Figure 32. Radiation mode coupling between two parallel interferometers: (a) step index 
profile; (b) BPM analysis. 
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Figure 33. Cross-sectional view at the center of LS (V = 0). 
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obviously impractical. The BPM provides an excellent method for analyzing the effects 
of the radiation mode coupling and is a useful tool for minimizing radiation mode effects 
while optimizing the interferometer spacing. The dependence of radiation mode coupling 
on interferometer spacing d, are clearly demonstrated through Figures 34(a) and 34(b). 
The interferometer spacing in the index profile of Figure 34(a) is reduced to d, = 5 um 
and the effects are directly demonstrated in the field distribution of Figure 34(b). 
Although extending the interferometers to infinity is not practical, it is anticipated that 
radiation coupling effects will drop off rapidly as d, is increased. The results of 
increasing the separation distance to d, = 20 um are clearly shown in Figures 35(a) and 
35(b). It is apparent from this representation that radiation mode coupling can be 
minimized utilizing the BPM as a design tool. 

1. Input/Output Characteristics 

The effects of radiation mode coupling can be quantified using the BPM 
algorithm. The two interferometers are initially spaced very close together. For each 
separation distance d, the input power to each interferometer and their respective output 
powers are recorded with V = 0. Due to the symmetry of the input power distribution 
equal power is delivered to the input of each interferometer. Figure 36 shows the input 
power and the output power for the left interferometer as a function of d, (normalized to 
the input power at the smallest separation distance). At small separation distances the 
output power is severely degraded due to the radiation coupling between the 


interferometers. The effects of the radiation mode coupling go through regions of 
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(b) 
Figure 34. Parallel interferometers with spacing d; = 5 um: (a) step index profile; 
(b) BPM analysis. 


(b) 
Figure 35. Parallel interferometers with spacing d; = 20m: (a) step index profile; 
(b) BPM analysis. 
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Figure 36. Input versus output power for the left interferometer. 
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constructive and destructive interference as d, is increased. However, these evaluations 
are performed with V = 0. The arm separation must be at least large enough to minimize 
the interference effects. 

2. Insertion Loss 

Insertion loss is also an important characteristic of optical interferometers. The 
induced AB due to the radiation mode interference tends to decrease the power 


transmitted and is reflected in the insertion loss calculation [5] 
Pr 
Li2(di) = 10 log sar fly (89) 
Pia 
where P, is the optical power that would be transmitted by the waveguide if the modulator 
were absent, and P,, is the intensity transmitted with the modulator in place and V = 0. 
To calculate the insertion loss a special index structure was constructed and is shown in 
Figure 37(a). The right interferometer has been replaced with a straight section of the 
waveguide. The magnitude of the optical field is shown in Figure 37(b). It can be noted 
from Figure 37(b) that even though the interferometer is replaced, the bends from the 
initial YPD and the branching angle cause perturbations in the optical field. The subtle 
insight obtained from this structure reinforces the fact that radiation losses have a direct 
dependence on branching angles. The insertion loss as a function of the separation 
distance d, is shown in Figure 38. As the distance between the interferometers increases, 
the insertion loss due to the induced AB decreases. Insertion loss of less than 1/2 dB is 


obtained at a separation distance of 4 um. For separation distances of greater than 4 um 


the radiation modes start to interfere constructively as was noted in the analysis of the 
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(b) 
Figure 37. Modified architecture for insertion loss and mode coupling calculations: 
(a) step index profile; (b) BPM analysis. 
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Figure 38. Insertion loss as a function of the separation distance. 
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output characteristics. This constructive interference increases the insertion loss 
somewhat and then tapers off at larger d,[5]. 
3. Radiation Mode Dependent Finite AB 
As previously mentioned, the induced AB due to the radiation mode interference 
gets smaller as d, gets larger. The expression that relates the transmission efficiency and 
the degradation in output intensity due to the radiation losses P,,,, is given as [5] 
Pioss = Pin cos?(AB L/2) — P12. (90) 
An approximation for AB can be found as [5] 
AB = (2) cos { Pit Pie) is (91) 
Since the input power changes for each d, due to the increase in path length, a 
normalization procedure is needed so that the induced AB can be monitored. First, the 
power loss through a stand-alone interferometer is computed, using the architecture of 
Figure 3 or Figure 37. The percent power loss incurred (assumed to be all radiation mode 
losses) is calculated to be p, = 5.8%. The output power at the given d; is then normalized 
using the respective input power and the loss value through the device p,. To normalize 
P,,,, to the input power at each separation distance we set 
Pioss = Pt Pin. (92) 
The approximation for AB is evaluated for several separation distances with the resulting 
AB shown in Figure 39. The induced AB drops off exponentially as the separation 


distance is increased. However the effects of the periodic interference are clearly shown 


in the 10 um to 14 um separation range. 
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Figure 39. Delta Beta as a function of the separation distance. 
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4. Modulation Depth 

The modulation depth of the individual interferometers and the corresponding 
effect that the radiation modes have on these modulation depths is also of concern. As 
mentioned earlier, the Input/Output analysis did not include the effects of radiation 
coupling with electrode voltages applied. The radiation-actuated AB reduces the 
modulation depth below 100%. A small amount of radiation interference can be tolerated 
if the performance of one interferometer does not adversely affect the other 
interferometer. Using the theoretical output intensity of Equation (63) for the 
symmetrical MZI the first electrode voltage point that produces output intensity 


extinction occurs at 


Aw, 


Vise 
* 2503303 


(93) 
where [ , r3;, and » are the same parameters defined for the single symmetric 
interferometer analyzed previously. The theoretical V, is calculated to be 8.17 V for the 
structural parameters given in Table 3. The BPM is used to examine the modulation 
depth of the left interferometer and the effect that the electrode voltage has on the 
adjacent interferometer as well. Figures 40, 41, 42, 43, and 44 show both the theoretical 
modulation characteristics and the BPM calculated outputs for the left and right 
interferometers for d; = 0.66 um, 1.33 um, 2.0 pm, 4 pm, and 10.66 um, respectively. 
Figures 40, 41, and 42 demonstrate the adverse effects of the radiation mode coupling 


from the left interferometer into the right interferometer. At greater separation distance 


the behavior closely follows the theoretically expected performance. However, the 
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Figure 40. Modulation characteristics (di + 0.66 um). 
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Figure 41. Modulation characteristics (di = 1.33 4m). 
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Figure 42. Modulation characteristics (di = 2.04.m). 
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Figure 43. Modulation characteristics (di = 4.0i.m). 
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Figure 44. Modulation characteristics (di = 10.66um). 
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modulation depth does not quite reach the zero limit predicted by theory. The output 
intensity of the right interferometer is also affected severely at small separation distance, 
even though no electrode voltage is applied to the right interferometer. At larger distance 
. the radiation mode coupling causes only a slight rise in the output intensity. 

The effects of the radiation mode coupling with applied electrode voltages are 
readily apparent when a cross section of the output optical field distribution is viewed. 
Figure 45 shows a cross-sectional view at the output, Z7 section, of the parallel array with 

4, applied to the left interferometer and V = 0 applied to the right interferometer. The 
mode power in the left interferometer has been coupled into the radiation modes as shown 
by the increased radiation modes adjacent to the left interferometer vice the radiation 
modes that exist due to the right interferometer. The effect of the absorption window is 
also readily apparent in Figure 45 as the radiation at the boundary is absorbed 
continuously. The total impact of the increased radiation modes due to applied electrode 
voltages is seen by comparing the view of Figure 45 to Figure 46, which is the same 


cross-sectional area with V = 0 on both interferometers. 
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Figure 45. Cross-sectional view at the output of the array (V = 8.17 volts). 
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Figure 46. Cross-sectional view at the output of the array (V 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


A. CONCLUSIONS 

The BPM has been previously utilized for integrated optical device analysis [4], 
however, the applicability of the 2-dimensional effective index method to these structures 
had not been shown empirically. The validation of the effective index approximation 
using the BPM through comparison to lab recorded data of a physical device was a major 
validation of this method. The BPM has been demonstrated as not only an effective 
analysis tool, but a potentially useful design tool as well. 

The development of the prefolded global propagator has made the BPM not only 
applicable, but computationally efficient. The BPM algorithm is written in a C™ code 
that can be ported to personal computer platforms. The memory and computational time 
requirements have been reduced sufficiently to make this analysis an attractive method of 
optical circuit analysis. 

B. RECOMMENDATIONS FOR FUTURE RESEARCH 

Although the BPM has been developed into a very effective simulation tool, the 
effective index method does not take into account variations in the vertical transverse 
plane that may have significant effects. New technologies in integrated optics are now 
taking advantage of laser trimming technology to modify slightly the physical structures 
of optical devices in order to make fine adjustments to inherent phase bias during the 
post-fabrication phase [10]. The asymmetric device modeled in this thesis was phase 


tuned by laser ablation as shown in Figure 47. The laser trimming is performed after 
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Ablation area 


Figure 47. Asymmetric interferometer showing laser ablation region. 
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manufacture in order to optimize linearity [10]. Due to variations in temperature or other 
imprecise fabrication results the inherent phase bias desired in such a device may not be 
exact. therefore some method of post manufacture processing is highly beneficial. The 
-BPM could easily be extended to a 3-dimensional structure in order to model these types 
of devices. The global propagator scheme would still be applicable to the 3-dimensional 
scheme so that computational requirements would only increase by the number of cross 
sections in the vertical grid. The power distribution through 3-dimensional devices 
would be viewed through cross-sections or the vertical dependence could be integrated 
out after the analysis so that total optical distributions could be viewed, such as the views 
generated in this thesis. 

The BPM could also be utilized to determine the feasibility of electrooptic 
switches that take advantage of the linear Pockels effect. A proposed design of a digital 
electrooptic switch is shown in Figure 48. A switching device such as Figure 48 could 
easily be simulated using the BPM method, and may have potential uses in future 
integrated optical components. Many potential applications of optical devices are 
constantly emerging and the BPM may have many uses in design and analysis of these 
devices. The key to the success of the BPM analysis on these new technologies are 


validity, speed of computation, and flexibility. 
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Figure 48. Proposed digital electrooptic switch. 
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APPENDIX. COMMAND LINE SCRIPT FILES 


This appendix contains some examples of typical command line script files that 
were used to run extended iterations of BPM analysis on the Sun workstations. 
1. EXAMPLE FOR SINGLE SYMMETRIC CODE 
This file passes lengths for L2 and L3 which are in microns. The 410 is for 
section L4, the 120 is for L4P coupler length, the 1000 is the electrode or L5 length, the 0 
0 are the left and right voltages, and the last digit tells the code that the entire index and 
optical field data sets should be saved to disk. However, note that in this case although 
the right voltage is passed, it is not used, it was left in for ease of use when changing from 


single to parallel systems. 


Zsing.h 
#! /bin/csh -f 
bpmP 12 12 410 120 10000 0 1 


Table 4 below gives a breakdown of the input line items. 
12 | 13 | L4 


VoltsL| VoltsR | Save | 
lbpmP 12 | 12 | 410 0 0 1 


Table 4. Input items. 


L4P | L5 
120 | 1,000 


2. EXAMPLE FOR ASYMMETRIC CODE 

This is a much simpler case. Although Dr. Bulmer of the Naval Reasearch 
Laboratory gave me the arm separation of 80 microns [10], I built the code to take any 
arm separation distance and calculate what the L2 and £3 lengths would be to support that 


structure. So the parameters passed in for the example are: arm separation is 80 um, V = 
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4.00 volts, and the data is to be kept. Caution should be exercised in this case. These 
data sets can run over 20 MB per run. It can shut the whole network down if caution is 


not exercised. 


Zasymmetric.h 
#! /bin/csh -f 
bpmP 80 400 1 


3. EXAMPLE FOR PARALLEL CODE 
This set up is in the exact same format as for the single symmetric case. The one 


item that should be observed in particular is the voltage input. The voltage is multiplied 


times 10°, so that 400 = 4.00 volts. 


Zparallel.h 
#! /bin/csh -f 
bpmP 1630 1000 1595 440 1000 001 


4. EXAMPLE FOR MULTIPLE RUNS 
The file name is left off the following example. The only real constraint in file 
names on the Unix system is that the file must have the executable flag enable in the file 
properties section. Typically "Z.h", "Zi.h" or some short variation was used. The 
following example makes 6 runs of the code and varies the left interferometer voltage 
from 0:0.1:0.5 volt. 
#! /bin/csh -f 
bpmP 1630 1000 1595 440 1000 0 00 
bpmP 1630 1000 1595 440 1000 10 00 
bpmP 1630 1000 1595 440 1000 20 090 
bpmP 1630 1000 1595 440 1000 30 00 
00 
00 


bpmP 1630 1000 1595 440 1000 40 
bpmP 1630 1000 1595 440 1000 50 
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